Ab-initio Dynamics of Rare Thermally Activated Reactions 
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We introduce a framework to investigate ab-initio the dynamics of rare thermally activated re- 
actions. The electronic degrees of freedom are described at the quantum-mechanical level in the 
Born-Oppenheimer approximation, while the nuclear degrees of freedom are coupled to a thermal 
bath, through a Langevin equation. This method is based on the path integral representation for 
the stochastic dynamics and yields the time evolution of both nuclear and electronic degrees of free- 
dom, along the most probable reaction pathways, without spending computational time to explore 
met ast able states. This approach is very efficient and allows to study thermally activated reactions 
which cannot be simulated using ab-initio molecular dynamics techniques. As a first illustrative 
application, we characterize the dominant pathway in the cyclobutene^butadiene reaction. 



Thermally activated conformational and chemical re- 
actions drive a large class of phenomena in physics, chem- 
istry and biology. From the theoretical standpoint, the 
understanding of such processes involves the characteri- 
zation of the ensemble of statistically significant reaction 
pathways, i.e. the set of consecutive system's configura- 
tions which are most likely to be visited during the tran- 
sition from the reactant to the product state, at a given 
finite temperature. Such pathways provide information 
about the structure of the transition state. 

In this context, the study of the dynamics of the reac- 
tion within an ab-initio approach is particularly challeng- 
ing. In fact, a formalism should be adopted which com- 
bines the quantum description of the electronic degrees 
of freedom with the out-of-equilibrium dynamics gener- 
ated by the coupling of the nuclear degrees of freedom 
with the heat bath. Ab-initio Molecular Dynamics (MD) 
simulations satisfy this requirement, at least in principle. 
In practice, all MD methods become very inefficient as 
soon as the activation barriers are of the order of few ksT 
or larger. The reason is that most of the computational 
time is spent to explore meta-stable states. 

In view of such difficulties, a variety of different meth- 
ods have been proposed, which neglect the dynamics aris- 
ing from the coupling with the heat-bath, but provide 
information about the structure of the potential energy 
surface in the transition region. In particular, methods 
based on the search of the minimum energy pathways 
have been proposed, such as the Nudged Elastic Band 
(NEB) [1]. Such approaches are based on the mini- 
mization of the total energy of a virtual chain of molec- 
ular configurations connecting the reactant and product 
states. They provide an estimate of the activation barri- 
ers, but do not account for the dynamics and thermody- 
namics of the reaction. For example, there is no guaran- 
tee that a NEB path is representative of the ensemble of 
statically significant reaction paths at a given tempera- 
ture. Methods such as the temperature-accelerated MD 



f2^ allow in principle to assess long-time dynamics, but 
they only apply to first-order reactions, are based on the 
approximations of the transition state theory, and are 
mostly efficient for small barrier problems. 

The Dominant Reaction Pathway approach (DRP)[3l 
SI El E] is a recently developed rigorous framework to 
identify the most probable transition paths in systems 
described by the over-damped Langevin dynamics. The 
DRP approach is particularly efficient when applied to 
study thermally activated transitions, since it does not 
waste computational time to explore metastable states. 
In addition, it directly yields the most probable paths 
without the need to perform a statistical analysis of a 
representative sample of reactive trajectories. It also 
offers a natural definition of the reaction coordinate, 
parametrized in terms of the total Euclidean distance 
covered by all atoms along the dominant trajectory. Fi- 
nally, the method yields information about the dynamics, 
since it provides the most probable time evolution of the 
system's degrees of freedom, during the reaction. 

So far, the DRP approach has been applied to study 
conformational reactions, in particular protein folding ^ 
6 , using empirical force fields. In this Letter, we extend 
this method and present its first application to an ab- 
initio simulation of a rare transition. This is done by eval- 
uating the molecular energy which enters the DRP equa- 
tions directly from quantum-mechanical calculations. As 
a result the DRP method yields the most probable nu- 
clear trajectories and the corresponding electron densi- 
ties, in the Born-Oppenheimer approximation. This im- 
provement opens the door to a systematic exploration of 
the dynamics of the electron density in reactions charac- 
terized by large activation energies. 

Let us begin by briefly reviewing the DRP ap- 
proach. The over-damped Langevin equation describ- 
ing a system of N atoms of nuclear coordinates X(t) = 
(xi(t), . . . , X7v(t)), in contact with a thermal-bath at 
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temperature T reads: 

ax _ D 

~dt ~ 



(1) 



where D is the diffusion coefficient, U (X) is the molecular 
potential energy, and r]{t) is a random noise with Gaus- 
sian distribution, zero average and correlation given by 

The conditional probability density for the system to 
be found in the configuration X/ at time t/, provided it 
was prepared in some initial configuration X^ at time ti is 
the Green's function of the Fokker-Planck operator asso- 
ciated to Eq. ([T]) and can be represented in the following 
path integral form: 



L/(X^)-C/(X^) 



where Seff is called the effective action, and is given by 



^eff 



[X] = dr^-^^Veff{X{r)). (3) 



Ve//(X) is called the effective potential, and reads 



K//(X) 



D 
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(4) 

If Xj and X^ are product and reactant configura- 
tions, respectively, then the factor exp(— S'e//[X]) inside 
the path integral expression ([2| represents the statistical 
weight of a given reactive trajectory X(t). In particular, 
the most probable trajectories are those which minimize 
the effective action functional Sef / [X] . These are the so- 
lutions of the classical equations of motion generated by 

the effective Lagrangian C = + ^eff with 

boundary condition X(ti) = X^ and X(tj) = Xj. 

The numerical advantage of the DRP approach 
follows from observing that such equations of mo- 
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tion conserve the effective energy E^ff 
Veff{^{t)). This property allows to switch from the 
time-dependent Newtonian description to the energy- 
dependent Hamilton- Jacobi (HJ) description. In the 
HJ framework the most probable pathways connect- 
ing the given initial and final configurations are ob- 
tained by minimizing the target HJ functional Shj[^] = 

/x/ dl yj^ [Eeff + Veff (X(0)], where dl = VdX^ is the 
elementary Euclidean distance in configuration space, 
along the dominant reaction path. 

The effective energy parameter E^ff fixes the total 
time ttot = tO^f) foi" a single barrier-crossing transition. 
In fact, the time t [X] at which the configuration X is vis- 
ited, during the most probable reaction pathway is given 
by: 



t 



(X)= / dl 



(5) 



distance C3 - C4 / A 



FIG. 1: Free-energy diagram for the reaction cyclobutene 
butadiene. Superimposed the NEB pathway (squares) and 
DRP pathway (circles), evaluated at T = 423K. 



In the present work, we have chosen E^ff = — ]/e//(Xj) 
which yields the longest total time ttot \^ LSj . We stress 
the fact that ttot is much shorter than the mean-first- 
passage time, as it corresponds to the time to reach the 
product, once the system has left the reactant. 

In practice, finding the dominant reaction pathway 
amounts to minimizing a discretized version of the ef- 
fective HJ functional: 



E 



[E^ff^Veff (X,)] A/,, 



where A^^ is the number of path discretization steps, and 
the effective potential Ve//(X) is determined according to 
Q from the molecular potential energy U (X) . The latter 
corresponds to the lowest eigenvalue of the Schrodinger 
equation for the electronic degrees of freedom, in which 
the nuclear coordinates are held fixed at the configuration 
X. Ali^i-^i is the Euclidean distance between the slices i 



and i + 1, i.e. A/^^^+i = y |Xi+i - X^ | . 

The global minimization of the HJ effective action (|6| 
is a delicate task. The difficulties arise from the rugged- 
ness of the effective potential, which depends both on 
the gradient and on the laplacian of the molecular po- 
tential /7(X). In order to perform such minimization, we 
tried different approaches, including simulated annealing 
and first-order methods such as conjugate gradients or 
the Broyden-Fletcher-Goldfarb-Shanno method [7]- The 
method which has performed best among the ones we 
tried is the Fast Inertial Relaxation Engine (FIRE) [8]. 
The noise in the numerical gradients of the HJ action be- 
comes a limiting factor, when the gradient norm descends 
below a certain threshold. Therefore, we employed a sim- 
plex algorithm for the final stage of the minimization. 

In the discretized representation of the HJ effective ac- 
tion (|6|, the width of the distribution of the Euclidean 
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FIG. 2: Reaction barrier evaluated along the dominant path- 
ways at different temperatures. The triangle represents the 
NEB result, within the same semi-empiric approach. 

distances between consecutive path slices, A/^^^+i, should 
not be allowed to increase in an uncontrolled way, in order 
to prevent all frames to collapse into the reactant or prod- 
uct configurations. A simple way to achieve this goal is to 
include a penalty function based on the standard devia- 
tion of the distances [9l, i.e. Sjjj [X] ^ 5'i/j[X] + C[X], 
with 

C[X(0] = A((A^2)-(A0')^A<72, (6) 

where (A?) = ^^Y.f=i^ ^ki+i, W = 
AT -1 Ylf=i^ and A is an adjustable pa- 
rameter. The penalty function limits the fluctuations of 
the Euclidean distances A/^^^+i and becomes irrelevant 
in the limit of very large number of path steps A^^. Un- 
fortunately, for a typical number of slices Ng 20 — 50, 
this term was seen to introduce significant bias. In order 
to overcome this difficulty, we introduced a Lagrange 
multiplier in the FIRE algorithm, which holds fixed at 
0.1 the ratio between the mean-square deviation from 
the average of the inter-slice distances and the square 
of the average inter-slice distance (A/). 

Let us now discuss the ab-initio evaluation of the 
molecular energy U (X) and of its first and second deriva- 
tives, which enter in the definition of the effective po- 
tential Q. This task requires to compute the ap- 
proximate ground-state energy of the Schrodinger equa- 
tion for the electronic degrees of freedom, in the Born- 
Oppenheimer approximation. In this first exploratory 
work, we choose to keep calculations as short as possi- 
ble and adopted the PM3 semi-empirical approach in the 
MOPAC implementation [10^. The pathway so obtained 
may be used as a starting point for more accurate calcu- 
lations, based, e.g., on density functional theory. 

As a first application of the ab-initio DRP method, 
we computed the dominant pathway for the reaction cy- 
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FIG. 3: Relative molecular energy along the DRP (dashed) 
and NEB (solid) paths, parametrized in terms of the curvilin- 
ear abscissa /, which corresponds to the total distance covered 
by the atomic nuclei during the reaction. 

clobutene butadiene shown in Fig{l] at the tempera- 
ture T = 423 K. In the course of this reaction an energy 
barrier of ~ 2 eV is overcome [11 and therefore this pro- 
cess can hardly be studied with ab-initio MD techniques. 
The DRP method was implemented using 16 discretiza- 
tion slices and the identification of the dominant reaction 
path required 80 CPU hours on 3 GHz processors. We 
verified that the FIRE algorithm achieved convergence 
to the same minimum, irrespective of the starting path. 

In Fig. [l] we present the resulting dominant reaction 
path, projected on the plane defined by the distance be- 
tween C3 and C4 and the dihedral angle defined by CI, 
C2 and the two hydrogens bonded to them. The DRP 
path (circles) is compared with the NEB path (squares) 
and superimposed on the free-energy map computed from 
ab-initio meta-dynamics [12] simulations. 

First of all, we observe that the DRP at this tempera- 
ture and NEB paths are quite different. While the min- 
imum energy path cuts across the iso- free-energy lines, 
the most probable path takes a shorter route, remain- 
ing for a longer portion in the high free-energy region. 
This can be explained by the fact that the dominant 
paths are the result of a compromise between minimiz- 
ing the total length of the path distance and the factor 
^/E^fJ^^Veff\^^JJ)]^ which appears in the integrand. In 
particular, in the high temperature limit, the free energy 
landscape becomes flat and the most probable path re- 
duces to a straight line in configuration space. On the 
other hand, as the temperature is reduced, the V/7(X) 
contribution to the effective potential becomes more and 
more important, and the path should approach the NEB 
result. In order to further study the temperature de- 
pendence of the DRP paths we compare in Fig. |2] the 
difference between the highest molecular energy reached 
along the dominant path and the reactant molecular en- 
ergy, for different temperatures. We observe that this 
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FIG. 4: Relative time at which each configuration is visited. 
Super-imposed are the electron iso-density surface lines, on 
the plane selected by the C1-C3-C5 atoms. The absolute value 
of the total barrier-crossing time Uot depends on the diffusion 
coefficient, and therefore on the details of the experimental 
conditions in which the reaction is realized. 



quantity grows linearly with temperature. A simple lin- 
ear extrapolation to zero temperature yields an energy 
barrier of 1.72 eV, which is close to the NEB result calcu- 
lated within the same semi-empirical approach, 1.76 eV. 
These results show that thermal effects can sizeably af- 
fect the structure of the statistically important paths at 
finite temperature. In particular, the most probable path 
does not in general pass through the saddle point of the 
molecular energy. 

It is also instructive to compare the value of molecu- 
lar energies along the NEB and DRP paths, reported in 
Fig. [Sj The NEB curve displays a small secondary maxi- 
mum, which is usually interpreted as a rotational barrier 
for the unsaturated hydrogen [11]. On the other hand, 
our DRP result does not show this feature, in the statis- 
tically most relevant path. This implies that the opening 
of the ring and the rotation of the unsaturated hydrogen 
are most likely to occur simultaneously, given the small 
value of the rotational barrier (about 0.01 eV) obtained 
at the PM3 level. 

The ab-initio DRP approach does not only yield path- 
ways in configuration space, but also allows to determine 
the time evolution of both nuclear and electronic degrees 
of freedom. As an illustration. Fig. [4] shows the rela- 
tive time t(l)/ttot at which each configuration is visited. 
Super-imposed are the electron iso-density surface lines, 
on the plane selected by the C1-C3-C4 atoms, at four 
different instants of the reaction. 

In conclusion, in this work we have developed an ab- 
initio approach to investigate the dynamics of molecular 
systems in contact with a heat-bath, and to find the most 
probable reaction pathway. Its main advantage is that 
its computational cost does not depend on the height 
of the barrier. Thus, for thermally activated reactions. 



the resources required are drastically smaller than that of 
ab-initio MD algorithms such as Car-Parrinello. This im- 
provement opens the door to systematically explore the 
time evolution of the positions of the atomic nuclei and 
of the electron densities during reactions characterized 
by large activation energies, and can be used to study 
the temperature dependence of the reaction mechanism. 
As a first application, we have applied our method to 
investigate a the cyclobutene butadiene reaction. We 
have found that, that the coupling with the heat-bath 
can appreciably affect the structure of the statistically 
important paths, at finite temperature. The information 
which can be obtained from DRP is therefore comple- 
mentary to that accessible from NEB. In fact, the latter 
can accurately predict the location of the saddle-point, 
but cannot investigate the dynamics of the statistically 
significant transitions. 

Calculations were performed on the WIGLAF cluster 
of the Physics Department of Trento University. The 
DRP method was developed in collaboration with M. 
Sega and H. Or land. We acknowledge important dis- 
cussions with W.J. Lester. 
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